home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / fitlsq.i < prev    next >
Text File  |  1995-07-26  |  3KB  |  83 lines
  1. /*
  2.    FITLSQ.I
  3.    Least squares fit a piecewise linear function to data.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func fitlsq(y, x, xp, weight=, stats=)
  11. /* DOCUMENT yp= fitlsq(y, x, xp)
  12.             ...
  13.         yfit= interp(yp, xp, xfit)
  14.      performs a least squares fit to the data points (X, Y).  The input
  15.      XP are the abcissas of the piecewise linear function passing through
  16.      (XP, YP) which is the best fit to the data (X,Y) in a least squares
  17.      sense.  The XP must be strictly monotone, either increasing or
  18.      decreasing.  As for interp, the piecewise linear fit is taken to be
  19.      constant outside the limits of the XP.
  20.  
  21.      The result YP is linearly interpolated through any consecutive
  22.      intervals of XP which contain no data points X, and linearly
  23.      extrapolated beyond the extreme values of X (if any of the intervals
  24.      of XP lie outside these extremes).
  25.  
  26.      A WEIGHT keyword of the same length as X and Y may be supplied in
  27.      order to weight the various data points differently; a typical
  28.      WEIGHT function is 1/sigma^2 where sigma are the standard deviations
  29.      associated with the Y values.
  30.  
  31.    SEE ALSO: interp
  32.  */
  33. {
  34.   np1= numberof(xp)+1;
  35.  
  36.   /* bin the input data into the xp, and create an extended version
  37.      of xp for which the bins may be directly used as an index list */
  38.   l= digitize(x, xp);
  39.   dx= (xp(0)-xp(1))*1.e30;
  40.   xx= grow(xp(1)-dx, xp, xp(0)+dx);
  41.  
  42.   xl= xx(l);
  43.   xu= xx(l+1);
  44.   dx= xu-xl;
  45.   g= (x-xl)/dx;
  46.   h= (xu-x)/dx;
  47.  
  48.   if (is_void(weight)) weight= 1.0;
  49.   hy= histogram(l, h*y*weight, top=np1);
  50.   hh= histogram(l, h*h*weight, top=np1);
  51.   gh= histogram(l, g*h*weight, top=np1);
  52.   gg= histogram(l, g*g*weight, top=np1);
  53.   gy= histogram(l, g*y*weight, top=np1);
  54.  
  55.   diag= hh(2:0)+gg(1:-1);
  56.   rhs= hy(2:0)+gy(1:-1);
  57.  
  58.   /* the triadiagonal system will be singular if there are any pairs
  59.      of consecutive bins -- remove these first */
  60.   list= where(diag);
  61.   diag= diag(list);
  62.   rhs= rhs(list);
  63.   off= gh(list)(2:0);  /* sub and super diagonal */
  64.   xpp= double(xp(list));
  65.   yp= TDsolve(off, diag, off, rhs);
  66.  
  67.   /* special treatment if endpoints removed allows linear extrapolation
  68.      of the fit until the true endpoints of the xp */
  69.   if (list(1)>1) {
  70.     yp= grow(yp(1)+(xp(1)-xpp(1))*(yp(2)-yp(1))/(xpp(2)-xpp(1)), yp);
  71.     xpp= grow(xp(1), xpp);
  72.   }
  73.   if (list(0)<numberof(xp)) {
  74.     grow, yp, yp(0)+(xp(0)-xpp(0))*(yp(0)-yp(-1))/(xpp(0)-xpp(-1));
  75.     grow, xpp, xp(0)
  76.   }
  77.  
  78.   /* add back any removed points */
  79.   if (numberof(xpp)<numberof(xp)) yp= interp(yp,xpp, xp);
  80.  
  81.   return yp;
  82. }
  83.